Spatial data

Lots of data in our field(s) are inherently spatial

R has lots of tools for interacting with and analyzing spatial data

the sp package

The sp package is the workhorse of spatial mapping and analysis in R.

It defines spatial classes (e.g. SpatialPoints, SpatialPointsDataFrame) that make it possible to work with spatial data.

actual spatial data versus coordinates

We often just have coordinate data, but this works differently in R from fully spatial data

myCoords <- data.frame(long = runif(20, min=35, max=36), lat = runif(20, min = 3, max=5))
plot(myCoords)

making spatial objects

You can use the SpatialPoints() function to make your coordinates into a fully spatial object

Warning: This function will assume that your columns are in the order: X coordinate, Y Coordinate

For a gold star: Which is the x and which is the y when we are dealing with latitude and longitude?

library(sp)
spMyCoords <- SpatialPoints(coords = myCoords)

Now there are a plethora of spatial things you can do with these points that you couldn't before

bbox(spMyCoords)
##            min       max
## long 35.009496 35.923433
## lat   3.079992  4.984301

making spatial objects

plot(spMyCoords, axes=TRUE)

maps

The maps package provides basic maps that can serve as a backdrop to your points

library(maps)
world <- map(database = "world")

maps

This is not an accurate election map!

USA <- map(database="state", fill=TRUE, col=c("red","blue"))

maps

kenya <- map(database = "world", region="Kenya")

Challenge

Plot the spCoords points as solid blue filled circles on top of the Kenya map

RGoogleMaps

library("RgoogleMaps")
UNCG <- GetMap(center = c(36.0674908,-79.8098429), zoom = 17)
PlotOnStaticMap(UNCG)

RGoogleMaps

UNCGSatellite <- GetMap(center = c(36.0674908, -79.8098429), zoom = 17, maptype = "satellite")
PlotOnStaticMap(UNCGSatellite, lon=-77.0493, lat=38.899, pch=16, col="red", cex=2)

Challenge

Get a satellite image of Africa, and plot a blue + symbol on Khartoum

GIS in R

You can do almost anything in R that you can do in ArcGIS.

This is accomplished through additional packages, notably rgdal which is an R interface to the Geospatial Data Abstraction Library (GDAL), which is a great piece of open source software for spatial analysis. GDAL is the brains behind QGIS.

The raster package deals with geospatial raster data, relying heavily on GDAL, as well.

The spatstats package implements a variety of spatial statitics.

rgdal

rgdal can read in most any vector data

The readOGR() function reads in vector data

First, download the shapefile from http://uncg.rstats.ninja/datasets/spatial/africa_simple.shp.zip. Then uncompress the file.

library(rgdal)
africa <- readOGR("/Users/wabarr/Dropbox/UNCG-workshop/datasets/spatial/africa_simple.shp",
                  layer="africa_simple")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/wabarr/Dropbox/UNCG-workshop/datasets/spatial/africa_simple.shp", layer: "africa_simple"
## with 1 features
## It has 3 fields

rgdal

spatial data can be plotted with the base R plot() function

plot(africa, col='maroon')

raster package

projections

proj4string(landCover)
## Loading required package: raster
## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(africa)
## [1] "+proj=eck6 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

Notice that the africa shapefile is in a different projection (Eckert VI equal area) rather than latitude / longitude, which is the most common coordinate system you will see

reprojecting

You can reproject vector data with sp::spTransform()

You can reproject raster data with raster::projectRaster()

reprojecting

Let's reproject the africa shapefile into lat/lon coordinates

africaLatLon <- spTransform(africa, proj4string(landCover))

plot together

Now that they are in the same coordinate system we can plot them together

plot(landCover)
plot(africaLatLon, add=T, border="purple", lwd=5)

mask raster with shapefile

maskedLandCover <- raster::mask(landCover, africaLatLon)
plot(maskedLandCover)

load in some point data

sites <- read.table("http://uncg.rstats.ninja/datasets/spatial/africanSites.txt", 
                    header=T, sep="\t")
sitesSP <- SpatialPointsDataFrame(coords=sites[,1:2], data=sites["siteID"])
plot(maskedLandCover)
plot(sitesSP, add=T, pch=16, col="blue")

buffers

Make 1 degree buffers using the `rgeos:gBuffer

rgeos::gBuffer(sitesSP, byid = T, width=1)
## class       : SpatialPolygonsDataFrame 
## features    : 5 
## extent      : -5.509803, 27.91958, -22.9821, 27.97751  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## variables   : 1
## names       : siteID 
## min values  :  site1 
## max values  :  site5